program sdd
use gra_forward_2
character*300 meshfile,denfile,ndenfile
character*100 str
integer ne,nn,nv
integer north,south,east,west,top,bottom
integer nut
integer i,j,k
real e0,n0,v0
real,allocatable::de(:),dn(:),dv(:)
real,allocatable::den(:)
real,allocatable::den3(:,:,:),nden3(:,:,:)
nc=iargc()
if(nc.ge.9)then
call getarg(1,meshfile)
call getarg(2,denfile)
call getarg(3,ndenfile)
call getarg(4,str) ; read(str,*) north
call getarg(5,str) ; read(str,*) south
call getarg(6,str) ; read(str,*) east
call getarg(7,str) ; read(str,*) west
call getarg(8,str) ; read(str,*) top
call getarg(9,str) ; read(str,*) bottom
else
        write(*,*) 'error parameter'
endif
call read_mesh_size(meshfile,ne,nn,nv)
     allocate(de(ne),dn(nn),dv(nv))
     allocate(den(ne*nn*nv))
     allocate(den3(ne,nn,nv))
     allocate(nden3(ne+1,nn+1,nv+1))
     
call read_mesh_dat(meshfile,ne,nn,nv,e0,n0,v0,de,dn,dv)

call read_phys_dat(denfile,ne,nn,nv,den)
write(*,*) 'Read data end'
!where(abs(den).le.0.1) den=0.0

call v2m(ne,nn,nv,den,den3)

call den2dden(ne,nn,nv,den3,nden3,north,south,east,west,top,bottom)

write(*,*) 'Calc data end'

call open_or_not(nut)
open(nut,file=ndenfile)
do i=1,ne+1
  do j=1,nn+1
    do k=1,nv+1
      if(nden3(i,j,k).ne.0)then
        write(nut,*) e0+sum(de(1:i-1)),n0+sum(dn(1:j-1)),v0+sum(dv(1:k-1)),nden3(i,j,k)
      endif
    enddo
  enddo
enddo
close(nut)

write(*,*) 'Process data Finished'

deallocate(de,dn,dv)
deallocate(den,den3,nden3)
end program
